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I INTRODUCTION 

One  of  the  major  problems  with  phenomenology  codes  for 
nuclear  explosions  is  our  lack  of  understanding  of  the  interaction 
between  multiple  bursts.  To  study  multiple  bursts  in  which  two  fire- 
balls interact  requires  a two-dimensional  radiation-hydrodynamics 
code  for  the  radiative  phase;  in  addition,  it  is  most  helpful  if  the 
later  hydrodynamics  can  be  studied  by  the  same  code,  rather  than  by 
transferring  back  and  forth  between  two  codes.  Several  such  efforts 
are  now  in  progress.  Extending  the  MICE  hydrodynamic  code  to  treat 
radiation  is  an  attractive  method:  in  addition  to  the  advantage  of  a 

single  code  for  early  and  late  times,  distinctive  features  of  this  code 
make  it  both  fast  (because  the  difference  equations  are  Implicit)  and 
flexible  and  simple  to  use  (because  of  the  technique  of  time-step 
splitting).  Recognizing  this,  the  Defense  Nuclear  Agency  authorized 
Mission  Research  Corporation  to  extend  this  code  to  treat  radiation 
diffusion , 

The  basic  MICE  hydrodynamic  code  is  described  in  Ref.  1.  It 
is  a one-fluid,  Eulerian  code  that  works  in  any  of  several  geometries. 
The  treatment  of  magnetic  fields  in  the  basic  code  has  been  dropped  in 
this  version,  which  is  intended  only  for  low  altitudes.  This  version 
treats  radiative  energy  transfer  in  the  approximation  of  radiation 
diffusion  (or  local  thermodynamic  equilibrium)  - that  is,  in  the  approx- 
imation that  the  photon  mean  free  path  is  small  compared  to  distances 
over  which  the  temperature  changes  appreciably.  This  approximation  is 
sufficient  for  early  times  for  bursts  below  40  km  altitude.  TJie 
treatment  of  radiation  currently  assumes  symmetry  about  the  vertical 
axis  and  uses  two-dimensional  cylindrical  (or  one-dimensional  spherical) 
coordinates. 


II  PHYSICAL  FORMULATION 


A.  DIFFERENTIAL  EQUATIONS 

We  denote  mass  density  by  p,  fluid  velocity  by  v,  specific 
internal  energy  (erg/gm)  of  the  fluid  by  I,  fluid  pressure  by  P, 
temperature  by  T,  the  local  acceleration  due  to  gravity  by  g = - ge  , 
and  radiative  energy  flux  density  integrated  over  all  frequencies  by 
F^ad  (erg/cm*-soc) . The  conditions  for  mass  balance,  momentum  balance, 
and  energy  balance  of  a volume  element  of  fluid  are 
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where  is  the  loss  rate  of  thermal  radiation  per  unit  mass,  which 
is  discussed  later  in  this  subsection.  Radiation  pressure  and  radiant 
energy  density  are  neglected. 


Radiation  Diffusion 


During  early  time  in  many  low-altitude  fireballs,  most  of 
the  radiant  energy  transfer  is  by  photons  whose  mean  free  paths  are 
much  smaller  than  the  fireball  radius.  The  fractional  temperature 
change  over  a few  mean  free  paths  is  small  for  all  frequencies  that 
are  significant  at  the  given  temperature.  In  such  a region,  the 
specific  intensity  is  close  to  that  for  thermal  equilibrium  at  that 
temperature,  and  we  can  treat  the  deviation  from  the  equilibrium  value 
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as  a small  quantity  (approximation  of  local  thermodynamic  equilibrium 
(LIE),  or  radiation  diffusion).  This  small  quantity,  which  characterizes 
both  the  departure  from  a Planck  spectral  distribution  and  the  departure 
from  isotropy,  is  characterized  by  the  radiative  energy  flux  density. 


The  frequency  spectrum  of  the  radiative  flux  is,  in  LTE,  a known 
function  of  the  temperature  and  absorption  coefficient  of  the  fluid. 


F^dV  = 


4tt 
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1 

(k"  + Kg)p 


3B^ 

“gq.  Vq  dv , 


where  k'’(p,T,v)  is  the  mass  absorption  coefficient  at  frequency  v, 
corrected  for  stimulated  emission,  ic^  is  the  mass  scattering  coeffi-. 
cient,  and  B^(v,T)  is  the  Planck  function;  the  flux  is  greater  at 
frequencies  at  which  it  experiences  less  absorption.  One  therefore 
, <^ftds  only  the  integrated  flux  density,  and  the  integration  over 
frequency  is  done  in  the  definition  of  the  Rosseland  mean  opacity, 
Xj^(p,T)  (w’'ich,  like  K^,  is  in  cm^/gm) . In  this,  the  photon  mean 
free  path,  1/(k^  + K^)p,  is  averaged  over  frequency  with  the  weighting 
function  9i.,  /9T,  which  emphasizes  high  photon  energies  around 
hv  4kT.  The  integrated  radiative  flux  density  is  then 


16oT^ 

3K,^p 


VT, 


(4) 


where  o is  the  Stefan-Boltzmann  constant.  Since  we  use  p and  I 
as  primary  thermodynamic  variables,  we  write  this  as 

(5) 

(6) 

is  the  radiative  flux  density  per  unit  temperature  gradient. 


^rad  = - 


where 


K(p,T)  = 16oTV3iCj^(P,T)p 
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The  cooler  part  of  the  burst  region  is  far  from  LTE.  Using 
the  diffusion  approximation  in  this  region  is  inaccurate  and  can  cause 
difficulties  in  obtaining  a numerical  solution.  Therefore,  we  have 
diminished  the  LTE  term  of  Eq.  (3),  - by  using  an  artificially 

high  Rosseland  mean  opacity  for  temperatures  kT  < 0.7  ev. 

Radiation  Loss  From  the  Fireball 

It  is  important  to  account  in  some  way,  however,  for  the 
radiant  energy  that  leaves  the  heated  region.  We  have  allowed  for 
such  non-LTE  radiative  energy  transfer  by  adding  to  Eq.  (3)  a radiative 
loss  term,  - which  is  explained  in  this  subsection.  This  radi- 

ative loss  term  transfers  radiant  energy  from  a volume  element  to 
infinity  (outside  the  heated  region) , whereas  the  LTE  term  transfers 
radiant  energy  only  between  adjacent  volume  elements. 

The  spontaneous  emission  per  unit  mass  per  unit  time  per 
unit  solid  angle  emitted  in  frequency  interval  dv  from  an  element 
of  fluid  is 

j^dv  = K"B^dv, 

where  again  ic''(p,T,v)  is  the  mass  absorption  coefficient  corrected 
for  stimulated  emission,  and  B^(v,T)  is  the  Planck  function.  To 
escape  completely  from  the  fireball,  the  emission  traveling  in  direc- 
tion Q traverses  a path  L(u)  from  the  radiating  element  to  outside 
the  fireball  (see  Fig.  1);  the  fraction  that  survives  this  traversal 
is  exp[-  T^(u)],  where 

T^(u)  = r k"(p,T,V)  p df 
‘^L(u) 

is  the  optical  depth  for  frequency  v along  the  path  L(u) . 
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Figure  1.  Axial  cross  section  of  interacting  fireballs,  showing 
geometry  used  to  calculate  radiation  loss.  The  power 
emitted  by  a fluid  element  into  solid  angle  dn  about 
a direction  u is  assumed  to  travel  along  the  path 
L(u).  Before  escaping  from  the  composite  fireball,  this 
power  is  attenuated  by  a factor  exp[-T(u)].  (The  dotted 
squares  are  the  cells  used  in  a representative  two-burst 
problem.  The  composite  fireball  shape  results  from 
independent  bursts  at  different  times  and  places.) 

Stimulated  emission  is  included  here  by  combining  it  with  the  spon- 
taneous emission  that  caused  it;  it  is  correctly  included  in  the 
fraction  exp[-  T^(u)]  by  the  correction  for  stimulated  emission  in 
k'  in  the  above  equation. 


We  consider  the  frequency  interval  0 < hv  < 4.13  ev  because 
it  includes  most  of  the  radiation  escaping  from  a low-altitude  fireball; 
this  interval  includes  the  infrared,  visible,  and  near  ultraviolet 
(wavelengths  > 3000  A)  and  is  referred  to  here  as  "thermal."  In 
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integrating  over  frequency,  it  is  con>fenient  to  use  the  Planck  (or 
"emission")  mean  opacity  for  this  frequency  range,  defined  by 


-‘I 


4 . 1 3 ev 


if"(p,T,V) 


B^(v,T) 


dv 


//■" 


3 ev 


B^(V,T) 


dv. 


(7) 

The  spontaneous  thermal  emission  per  unit  mass  per  unit  time  per 
unit  solid  angle  emitted  by  a fluid  element  is  therefore 


i4(P,T) 


•♦.IS  ev 

k"(P,T,v)  B^(v,T)  dv. 


= ><P4(P>1')  otVtt,  (8) 

where  the  integral  of  the  Planck  function  over  this  frequency  range 
is  written  as  fj^(T)  oT'Vtt.  The  function  f/^(T)  is  taken  from  Allen's 
table  (Ref,  2) , To  find  how  much  of  this  thermal  emission  escapes 
from  the  fireball,  we  make  the  approximation  that 

,1 3 ev 

exp[-  T^(u)]  dv=s  exp[-  t(u)], 

where 


T(u)  H j Kj,^(p,T)pde  (9) 

L(Q) 

is  the  (Planck)  mean  optical  depth  for  this  frequency  range  along  the 
path  L(u) . 'Hiis  approximation  of  using  a single  attenuation  factor 
over  the  whole  frequency  range  becomes  less  serious  as  the  fireball 
grows  optically  thin. 
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Integrating  over  all  directions  then  gives  the  loss  rate 
per  unit  mass  of  thermal  radiation  (erg/gm/sec)  from  a fluid  element, 

\ad  ^ U f 

where  dJ2  denotes  an  element  of  solid  angle  about  the  fluid  element. 
The  factor  is  outside  the  integral  because  spontaneous  emission 
is  isotropic. 

Eqs.  (5),  (8),  (9),  and  (10)  replace  and  in 

Eq.  (3).  Then,  for  problems  having  symmetry  about  the  vertical  axis, 
the  system  of  Eqs.  (1)  to  (3)  consists  of  four  component  integro- 
differential  equations  for  the  variables  p,  1,  v^,  and  v^. 

B.  EQUATIONS  OF  STATE  AND  OPACITY 

To  complete  the  problem,  we  must  specify  algebraic  or  tabular 
relations  for  these  properties  of  air  of  sea- level  composition:  fluid 

pressure  P(p,l),  temperature  T(p,I),  Rosseland  mean  opacity  iCj^(p,T), 
and  Planck  mean  opacity  for  the  frequency  range  0 < hv  < 4.13  ev, 

Kp4(P,T). 


Equations  of  state  for  air  are  needed  to  determine  fluid 
pressure  P(p,I)  and  temperature  T(p,I)  in  terms  of  density  p and 
specific  internal  energy  I.  Our  equation  of  state,  written  by  D. 
Sappenfield,  uses  data  from  Ililsenrath  and  others  (Refs.  3 and  4)  for 
1500°  < T<  5x10®  °K  (0.13  < kT  <:  430ev) , data  from  Gilmore  (Ref.  5) 
for  T = 1000  °K,  and  relations  for  an  ideal  diatomic  gas  with  no 
vibrational  excitation  for  T < 500  °K.  TJic  zero  of  internal  energy 
refers  to  gaseous  Nz  and  O2  at  T = 0°K. 
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Data  tables  are  used  to  do  two-dimensional  linear  interpo- 
lation with  log  p,  Zn  I,  tn  T,  and  P/pI  as  variables.  To  save 
computation,  we  use  this  routine  to  make  larger  data  tables  (10  den- 
sities by  75  interna]  energies)  in  which  we  do  two-dimensional  linear 
interpolation  with  p,  I,  T,  and  P/pI  as  variables.  From  this  inter- 
polation we  also  get  the  derivatives  like  (31/31)^  for  Eq,  (5);  these 
are  piecewise  constant.  Figs.  2 and  3 present  the  two  equation  of 
state  relations  in  a more  conventional  form,  showing  I and  P/pI  as 
functions  of  temperature  for  several  densities. 

The  Rosseland  mean  opacity  of  air  is  calculated  from  an 
analytic  relation  of  Erode  (Ref.  6);  this  is  plotted  in  Fig.  4.  The 
maxima  near  kT  » 6 ev  occur  because,  for  kT  < 2 ev,  few  photons  have 
enough  energy  to  ionize  air,  while  for  kT  > 10  ev,  ionization  has 
greatly  reduced  the  number  of  bound  electrons.  Figure  5 shows  the 
corresponding  photon  mean  free  path,  Xj^(p,T)  = l/<j^(p,T)p,  Brode's 
relation  gives  an  artificially  high  opacity  at  low  temperatures 
(kT  < 0.7  ev)  for  the  reason  mentioned  on  page  7.  Computing  time 
for  this  formula  was  reduced  by  60%  by  recoding. 

The  Planck  mean  opacity  of  air,  averaged  over  the  photon 
energy  range  0 < hv  < 4,13  ev,  is  taken  from  Sowle  and  others  (Ref.  7, 
pg,  10)  in  the  form  of  tables  of  Kp^(p,T).  Again,  Fig,  6 shows  this 
opacity,  and  Fig.  7 shows  the  corresponding  mean  free  path, 

^P4(P,T)  = l/iCp^(p,T)p.  This  opacity,  like  the  Rosseland  opacity 
above,  includes  the  correction  for  .stimulated  emission.  As  before, 
these  data  tables  are  expanded  by  two-dimensional  linear  interpolation 
with  £n  p,  Zn  T,  and  £n  Kp^  as  variables;  this  larger  table  is  used 
to  do  two-dimensional  linear  interpolation  with  p,  T,  and  Kp^  as 
yariables , 


Density,  p (gm/cm^) 


I 


ROSSELAND  MEAN  OPACITY 


PLANCK  MEAN  OPACITY  FOR  0 to  4.13  ev  PHOTONS,  k^a  (cm^/gm) 


TEMPERATURE,  kT  (ev) 


Contour  plot  in  the  density-temperature  plane  of  the  mean 
free  path  of  photons  in  air  (based  on  the  Planck  mean 
opacity  for  photon  energy  range  0 to  4.13  ev).  This  plot 
shows  the  data  of  Figure  6 in  the  form 
^p^CPjI")  ” 1 P >T ) p. 


INITIAL  AND  BOUNDARY  CONDITIONS 


Undisturbed  Atmosphere 


At  low  altitudes,  the  dependence  of  mass  density  p(z) 
on  altitude  z is  from  the  1962  U-S.  Standard  Atmosphere  (Ref.  8). 

The  pressure  just  above  the  top  of  the  grid  is  also  taken  from  this 
source,  and  pressui’es  P(zj  at  lower  altitudes  are  obtained  by  inte- 
grating the  hydrostatic  equilibrium  equation.  The  acceleration  due  to 
gravity  is  g(z)  = goRo/(Ro  where  go  = 980.665  cm/sec^  and 

Ro  = 6.37817  X 10®  cm.  The  specific  internal  energy  is  then 
I(z)  = P(z)/0.4  p(z).  The  curvature  of  the  earth  is  neglected. 


Deposition  of  Initial  X-rays 


The  initial  X-rays  are  characterized  by  an  energy  spectrum 
j de) , defined  so  that  the  X-ray  energy  in  the  photon  energy  interval 
de  is  J(e)  de.  We  usually  use  a Planck  spectrum  characterized  by 
some  temperature  T and  total  energy  E , 

X X 


C e®  de 

= exp  (e/kT^)  -1 


for  e > 0.013  kev. 


where  the  constant  C normalizes  the  total  X-ray  energy: 


0.0 1 3 kev 


J(e)  de  = E^. 


The  X-rays,  assumed  to  come  from  a point  source,  could  be  further 
characterized  by  an  angular  distribution,  which  we  currently  take  to 
be  isotropic. 
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We  consider  energy  loss  of  these  X-rays  only  by  photon 
absorption]  Compton  scattering  is  ignored.  The  problem  is  to  calculate 
the  change  of  specific  internal  energy  AI  of  air  caused  by  deposition 
of  X-ray  energy  along  rays  from  the  burst  point  on  which  the  density  is 
given.  Our  deposition  routine,  written  by  D.  Glenn,  uses  the  method  of 
Sappenfield  and  Tierney  (Ref.  7,  pp.  3 to  8);  the  part  relevant  at  low 
altitude  is  described  b,iow. 

For  air  that  is  not  fully  ionized,  our  corrected  mass 
absorption  coefficient  is 


K"(e)  = 


for  e > 0.4  kev, 

for  0.013  < e < 0.4  kev. 


(13) 


where  ki  = 3620  cm^-kcv^/gm  and  ka  ■ 369  cm^-kev^'*^  /gm;  Fig,  8 
shows  this  absorption  coefficient.  The  first  part  of  this  relation 
is  based  on  data  of  McMaster  (Ref.  9);  the  second  part  is  close  to  data 
for  molecular  oxygen  in  the  DNA  Reaction  Rate  Handbook  (Ref.  10) . This 
rough  absorption  coefficient  includes  no  dependence  on  the  density  of  the 
air  and  none  (except  as  in  the  next  paragraph)  on  its  temperature. 


For  air  that  is  fully  ionized  by  X-rays,  we  assume  that  the 
resulting  photoelectrons  thermalize  in  a time  shorter  than  that  (a  fe\' 
microseconds)  during  which  the  initial  X-rays  are  emitted  by  the 
debris.  That  is,  the  photoelectrons  resulting  from  the  earliest  X-rays 
from  the  debris  quickly  undergo  collisions  that  free  other  electrons 
and  spread  the  energy,  returning  the  air  toward  thermal  equilibrium. 

When  all  electrons  in  a region  have  been  so  freed  and  so  no  further 
(net)  X-ray  absorption  can  occur,  the  specific  internal  energy  of  the 
air  there  has  reached  an  upper  limit.  This  limit  is  the  internal  energy 
at  which  air  close  to  thermal  equilibrium  is  nearly  fully  ionized,  or 
about  I = 3 X 10^^  erg/gm.  We  use  this  limit  in  determining  an 
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MASS  ABSORPTION  COEFFICIENT,  K^Ce)  (cm=^/gm) 


WAVELENGTH  X (A) 


Figure  8.  Mass  absorption  coefficient  of  incompletely  ionized  air 
vs.  photon  energy.  The  absorption  edge  is  at  the  energy 
of  the  K absorption  edge  of  nitrogen. 

effective  mass  absorption  coefficient  for  air  that  becomes  fully 
ionized: 

Tk  /e^  for  e > 0.4  kev, 

K*(e)  =)  (14) 

(o  for  e < 0.4  kev, 

* * 

where  k is  a constant  (k  < ki)  to  be  determined  by  the  condition 
that  the  specific  internal  energy  not  exceed  I (k  and  R will  be 
determined  together  through  Iteration  with  Eqs.  18  and  15).  For 
simplicity  we  assume  that  photons  of  energy  less  than  0.4  kev  are  not 
absorbed  in  the  fully  ionized  region;  these  photons  carry  only  a 
small  part  of  the  energy  from  most  sources  of  interest. 


We  determine  first  the  size  of  the  region  in  which  air 
becomes  fully  ionized.  For  simplicity  we  assume  that  this  region  is 
a sphere;  that  the  density  in  it  is  a constant,  the  current  burst- 
point  air  density  po ; and  that  the  specific  internal  energy  in  it 
initially  is  negligible  compared  to  I*.  The  radius  R*  of  this 
region  is  determined  by  balancing  the  energy  absorbed  in  this  sphere 
with  its  net  energy  intake: 


* 

I Po 


■I 


0.**  kev 


J(e)  [1  - exp  (-  K (e)  po  R*)]  de. 


By  assumption  (Eq.  14),  photon  energies  below  0.4  kev  do  not  contri- 
bute. 


Outside  the  fully  ionized  region,  we  calculate  the  X-ray 
energy  deposition  in  spherical  coordinates  (r,  6)  centered  at  the 
burst  point  and  later  transform  the  results  to  the  desired  grid.  The 
angle  0 is  measured  from  the  vertical.  In  this  region  the  unabsorbed 
X-ray  energy  in  energy  interval  de  that  crosses  unit  area  at  a 
distance  r is 


F^(r,0,e)  de  = exp 

4irr^ 


K (e)  po  R - P(t',0)  dr" 


0 


for  r > R , 


where  p(r,  0)  is  the  current  mass  density,  which  may  have  been 
altered  by  a previous  burst.  The  energy  per  unit  mass  deposited 


at  position  (r,  0)  is  then 

* 
I 


for  r < R , 


J K"(e)  Fg.(r,  0,  e)  de 


for  r > R . 


(1/ 


0.013  kev 


AI  (r,  0)  = 


Continuity  of  AI(r,  0)  at  the  boundary  requires  that 


•'A 


0.4  kev 


K-(e)  F^(R  , e,  e)  de; 


(18) 


the  contribution  to  this  integral  from  energies  below  0.4  kev  is  ignored 

★ ■* 

for  simplicity.  We  determine  k and  R simultaneously  from  conditions 
(18)  and  (15)  by  iteration;  after  that  the  deposition  of  X-ray  energy 
by  Eq.  (17)  is  straightforward. 

Our  deposition  outside  the  fully  ionized  region  uses  25  angu- 
lar zones  (A0  = 7°. 5)  and  50  nonuniform  radial  zones.  The  transforma- 
tion back  to  the  original  grid  uses  two-dimensional  linear  interpolation 

/r 

p(r'',0)dr 

is  evaluated  by  the  trapezoidal  rule.  k" 


Initial  Debris  Energy 


The  debris  is  characterized  by  its  mass  (if  this  is 
important)  and  its  initial  kinetic  energy  K^,  Large  velocity  grad- 
ients near  the  burst  point  cause  difficulties  in  our  numerical  method, 
however,  so  we  elect  to  convert  this  kinetic  energy  immediately  to 
internal  energy.  This  procedure  seems  acceptable  because  the  results 
show  an  appropriate  partition  of  the  total  energy  in  the  problem  between 
kinetic  and  internal  by  2 or  3 msec  after  the  burst.  The  debris  mass 
(if  any)  and  internal  energy  are  initially  deposited  uniformly  through 
a sphere  of  radius  Ro,  which  we  choose;  typically  Ro  is  a few  times 
the  initial  cell  size. 
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Boundary  Conditions 


As  long  as  the  outer  boundary  is  well  outside  the  fireball, 
radiation  diffusion  is  unimportant  there.  We  have  made  the  simple 
approximation,  which  has  been  satisfactory,  that  radiation  diffusion 
causes  , o heating  or  cooling  of  the  fluid  at  the  outer  boundary, 

V*Frad  = The  general  regridding  capability  enables  this  outer 
boundary  (in  the  shape  of  a cylinder,  or  sphere)  to  be  moved  to  keep 
it  well  outside  the  fireball (s). 

The  hydrodynamic  part  of  the  boundary  conditions  must  allow 
shocks  to  pass  through  the  outer  boundary  smoothly,  without  creating 
artificial  reflected  waves  at  the  boundary.  At  this  boundary  we 
assume  (1)  that  there  is  no  inflow  of  mass,  and  (2)  that  the  time 
rate  of  change  of  the  amplitude  of  certain  incoming  waves  is  zero  (see 
Ref.  1,  Appendix  C) . We  allow  shocks  to  pass  through  the  boundary 
when  they  become  moderate  to  weak,  and  these  vonditions  seem  to  give 
satisfactory  results . 
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III  MATHEMATICAL  TECHNIQUE 

A.  DIFFERENCE  EQUATIONS 

For  problems  having  symmetry  about  the  vertical  axis,  the 
system  of  equations  (1)  to  (3)  consists  of  four  component  integro- 
differential  equations  for  the  variables  p,  I,  v , and  v . Each  is 
first  order  in  time  and  space,  except  that  radiation  diffusion  makes 
Eq.  (3)  second  order  in  the  spatial  variables.  The  MICE  code  solves 
this  system  by  a modified  Lax-Wendroff  method;  the  method  uses  the 
technique  of  time-step  splitting  and  is  modified  to  be  implicit. 

Our  treatment  of  radiation  involves  further  time-step  splitting  (see 
Eqs.  19)  and  is  also  implicit  (except  for  the  radiation  loss  term, 
which  involves  integrals) . 

The  difference  equations  are  written  in  two-dimensional 
cylindrical  coordinates  (r,z) . A (nonuniform)  grid  of  J by  J 
spatial  cells  is  specified,  with  boundaries  rj^  ^^|k=l , . . . , J^+1)  and 
z.  (Jc=l,.,.,J  +1).  The  coordinates  of  the  cell  centers  are  denoted 
by  rj^  (k=l,...,J^)  and  (2.=!, . . . ,J^)  . The  choice  of  cell  boundaries 

is  arbitrary  (Sec.  III-C  gives  our  criteria);  however,  if  adjacent  cells 
differ  in  width  by  more  than  10  or  20  percent  for  a time  long  enough 
for  their  contents  to  interact,  the  results  are  degraded.  The  code  is 
limited  to  < 98  cells  vertically  (which  give  a maximum  of  about 

5000  ceils).  For  the  case  plotted  in  Fig.  5 and  for  those  mentioned  in 
Sec.  IV,  the  cell  size  in  the  firebail(s)  ranges  from  3x3  meter  at 
burst  time  to  20  x 20  mete,  (at  most)  after  a few  hundredths  of  a second. 


We  wish  to  make  ri  = 0,  so  we  treat  the  axis  cells  specially, 
defining 
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Tlie  dependent  variables  are  advanced  in  time  by  application 
of  four  difference  operators,  R , R , L , and  L , that  are  one-dimen- 
sional  (that  is,  each  operator  involves  differences  in  only  one  dir- 
ection, specified  by  its  subscript).  Tliese  operators  are  applied 
sequentially  instead  of  simultaneously  (time-step  splitting) . In 
describing  this,  we  denote  dependent  variables  such  as  p(rj^,z^,t*^)  by 
p|^  ^ etc.,  and  consider  the  collection  of  dependent  variables  at  all 
grid  points  at  time  t*^  to  be  a vector,  having  4 J J components: 

X z 


u"  = (p 


^Vk,r  fVk,£’ 


k— 1 , . • . >*^^5  "^"1  > » • • * 


We  define  three  successive  sets  of  intermediate  variables,  each  result- 
ing from  the  application  of  one  more  matrix  difference  operator: 


yn,r^  + R (u",y”’^). 


^n,r  ^ ^^^^n,r^yn,rz^^ 
y"-rzp=  un,rz  ^ j.yn,rz^^n,rzp^  _ 

Finally  the  last  operator  yields  the  time-advanced  variables: 


(19) 


„n,rzp  ^ . 

Without  the  latter  of  the  two  operands  in  parentheses,  each  such 
equation  would  represent  a collection  of  4J^J^  explicit  equations; 
with  it,  each  equation  is  implicit  and  represents  a system  of  4J  J 

X z 

simultaneous  algebraic  equations  for  that  same  number  of  variables. 

Basically,  changes  in  the  variables  due  to  hydrodynamic  terms 

on  the  right  sides  of  Eqs.  (1)  to  (3)  are  accomplished  by  the  and 

L operators,  and  changes  due  to  the  radiation  terms  of  Eq.  (3)  are 
z 


accomplished  by  and  R^.  and  are  the  "L^  phase  2"  and 
"L  phase  2"  operators  of  Fajen  (Ref.  1,  pp.  18-23);  his  "L  phase  1" 
operators  are  unnecessary  in  the  absence  of  a magnetic  field.  The 
R^  and  R^  operators  affect  only  the  j,  and  not  densities  or  veloc- 
ities (because  radiation  terms  appear  only  in  the  energy  equation, 

(3)).  Thus,  the  R difference  operator  consists  of  linear  difference 

^ n r n r 

equations,  one  for  each  cell  that  relate  £,  and 

^k+i  t known  quantities. 


The  R and 
r 


R^  operators  are  defined  together  (in  Eqs.  20 


and  21),  using  the  following  abbreviated  notation.  The  coordinate  x 

can  be  either  r or  z,  and  we  retain  only  the  subscript  denoting  the 

X coordinate.  For  example,  and  denote  ^ and 

if  X = 4;  if  X = z they  denote  I, ’.  and  I. ' . ' (the  second  super- 

K > 3 K , 3 

script  is  abbreviated  in  this  case).  The  cell  widths  and  the  distances 

between  cell  centers  are  denoted  by  Ax.  = x.  - x.  and 

} 3+1/2  3-1/2 

+ + respectively,  and  is  a .space  interpolation 

coefficient,  d.^,  = (x.^  , -x.)/(x.  ,-x.).  The  index  m used  in 
3 + ‘ 3 + 1/2  3 3 + ‘3 

representing  divergences  is  1 for  the  r direction  and  0 for  the 


z-direction. ' 


TJie  specific  heat  capacity  is  denoted  by  c^  = (SI/OT)^, 


and  functions  of  thermodynamic  variables  are  abbreviated  as 
Tj  = T(Ij,p^),  Kj  = K(Pj,T^),  etc.  A time  interpolation  coefficient 
0 allows  the  R operators  to  be  explicit  (0=0)  or  implicit  (O<0<1) ; 
we  use  0=1.  (We  tried  0=0  but  had  difficulties.) 


In  representing  the  radiative  flux  density,  note  that  in 
Eq.  (5)  the  second  term  in  the  brackets  is  generally  small  compared  to 
the  first  because  T(I,p)  depends  only  weakly  on  p (see  Fig.  2).  We 
therefore  treat  9p/9x  less  elaborately  than  91/9x  here.  The  parallel 
component  of  the  radiative  flux  density  F at  a cell  boundary  is 
represented  as 


In  spherical  coordinates  m - 2 for  the  r direction;  in  Cartesian  1 
coordinates  ra  = 0 for  all  directions. 
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I 
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and  the  coefficients  composed  of  quantities  at  cell  boundaries  ' 
are  evaluated  by  linear  spatial  interpolation, 


The  radiation  terms  in  the  energy  balance  (Eq,  3)  are 
complicated  by  the  integral  in  (Eq.lO).  This  is  treated  in 


the  next  subsection;  the  part  of  that  is  included  in  the 


operator  is  denoted  by  (D  ).,  given  by  Eq.(29).  Hie  effect  of 

X ■} 


the  radiation  terms  on  the  energy  balance  (Eq.3)  ic  then  repre- 
sented by 


in,x_in 

j ~ 3 


1 


x"*  , (F 

l-l/z  X'^l-i/2 


- fD  V? 


I 


At  the 
causes 
(k=J^) 


Two  boundary  conditions  complete  the  system  of  equations  (21), 
outer  boundary,  well  outside  the  fireball,  radiation  diffusion 

“V 

no  heating  or  cooling  of  the  fluid,  “ 0.  At  the  outside 

and  top  this  condition  is 


in,x  _ jn 


At""V2 


- (Vj 


(22) 


The  other  condition,  which  applies  at  the  bottom  (5.»1)  and  at  the  axis 
(k*l),  includes  the  symmetry  requirement  at  the  axis: 


jn,x  ,n 
M X 

At"^V2 


Pi  (T3/2  -0) 


(F 


3/2 


if  x = r, 
if  X = 2. 


(23) 


Equations  (21)  to  (23),  with  the  substitution  (20),  form  a set  of 
linear  algebraic  equations  in  the  unknowns  I^’  . The 

system  is  tridiagonal  and  is  solved  by  a simple  elimination  scheme. 
As  stated  above,  the  R operator  makes  no  change  in  densities  or 
velocities . 


B . EQUATION  FOR  .RADIATION  LOSS 

In  calculating  the  loss  rate  per  unit  mass  of  thermal 

radiation  (Eq.  10),  we  evaluate  the  angular  integral  by  a sum  over 
four  directions: 

Drad  ^ T4(P,T)  j^cxp(-x(r))  + exp  (-x(-r))  + exp  (-t(z))  + exp  (-T(-z)jj  AQ, 
where  AQ  = v sterad.  This  form  makes  it  convenient  to  split 
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and  include  the  first  two  terms  in  the  operator  and  the  last  two 
in  the  operator.  We  denote  the  two  parts  by  and  so 

that 


D , = D + D . 
rad  r z 


(24) 


We  can  then  use  the  abbreviated  (single-subscript)  notation  of  the 
previous  subsection,  in  which  x denotes  i or  z and  only  the 
subscript  referring  to  the  x direction  is  retained. 

In  evaluating  t(x)  for  a given  cell,  wc  use  Eq,  (9)  and 
assume  that  p and  I are  constant  throughout  each  cell.  The  optical 
thickness  (for  0 to  4.13  ev  photons)  of  cell  j along  a direction  x is 


At^(X)  . Kp^(pJ,T5)  Ax. 


for  1 < j < 


(25) 


j 


and  the  optical  depth  from  the  right  (or  top)  edge  of  cell  j to  the 
right  (or  top)  boundary  of  the  problem  in  the  x direction  is 


t"  , (x)  -V'  At"(x) 
J + l/2  Z-/  1 

i=j  + i 


for  0 < j < J -1,  (26) 

A 


and  that  from  the  left  (or  bottom)  edge  of  cell  j to  the  left  (or 
bottom)  boundary  of  the  problem  in  the  -x  direction  is 

J 


z * t 

i=J 


(?) 


T.  (-X)  - 
3-1/2 


1^1 

'z 

^i=i 


if  X = r. 


if  X = z. 


(27) 
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Because  of  the  coarseness  of  our  grid,  some  cells  that 
radiate  strongly  are  optically  thick  (At?  > 1) ..  In  order  to  account 
reasonably  for  the  self-absorption  within  a cell  we  allow  t(x)  and 
t(-x)  to  vary  through  the  cell. 


t(x)  = k„.|p?,T?)  p?  fx.  -x]  + T?  (x) 

j ^ 3+1/2  I 3 + 1/2  '•  ■' 


and  similarly  for  t(-x),  and  average  over  the  mass  of  each  cell. 

Thus,  the  part  of  that  is  included  in  the  operation  is 


(D  )?  = -^ 

x^'3 


Xe^e  exp[-T(-$c:jjA 


p dV 


f P dV 
*'ce££  3’ 


here  Cj^  is  a calibration  constant  discussed  in  Sec.  III-D.  We  assume 
that  p and  j^(p,T)  are  constant  throughout  the  cell,  and  we  ignore 
any  variation  through  the  cell  of  its  crossrscctional  area  perpendicular 
to  the  X direction.  An  integral  useful  in  evaluating  the  above 
equation  is 


jl-expr-AT?(x)  I 

ij-  J exp  [-TCS)]  dx  = exp[-T.^_^^ 


D-j/2 


At. ‘(5c) 
1 


and  the  final  result  for  (D  ) . is 

X’' j 


f 


I 


I 


I 


I > 

i 


L 


A simpler  calculation  (setting  x=x^  in  Eq,  28)  would  have  given 
exp[-0.5ATj (x)]  in  place  of  the  expression  in  large  braces;  for 
optically  thin  cells  (ATj(x)«l)  these  expressions  are  equal,  but 
when  one  cell  is  optically  thick  and  those  in  front  of  it  are  optically 
thin,  only  Eq.(29)  reduces  correctly  to  blackbody  radiation  from  the 
surface  of  the  cell. 

C.  CHOICE  OF  TIME  STEP  AND  SPATIAL  GRID 


An  implicit  difference  scheme,  such  as  used  here  to  treat 
radiation,  is  numerically  stable  with  much  larger  time  steps  than  an 
explicit  scheme.  The  present  scheme  for  radiation  diffusion  ii  stable 
enough  so  that  accuracy,  rather  than  numerical  stability,  limits  the 
size  of  time  step.  We  limit  the  time  step  by  requiring  that  the 
change  in  specific  internal  energy  of  each  cell  during  each  operation 
(R^  and  R^)  must  not  exceed  a given  fraction  c^^^  (we  use  c^^^  = 0.025) 
of  the  maximum  specific  internal  energy  anyAvhere  at  that  time: 


n+i/2  _ 
rad 


(30) 


At  later  times  numerical  stability  of  the  hydrodynamics  also  limits  the 
time  step  (Ref.  1,  pg.  2): 


with 


We  find 


Atr"''/"  = 


Ax. 


k,l 


At""V2  =minfAt""y^  At^V^ 

\ rad  hyd 


c.  , = 0.2  satisfactorv. 
hyd 


(31) 


(32) 
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The  spatial  grid  is  chosen  to  have  a finely  gridded  region 
of  21  X 41  cells  around  each  burst  during  its  growth  stage.  (When  a 
second  burst  occurs,  the  finely  gridded  region  is  achieved  by  having 
a vertical  strip  of  narrow  cells  intersect  a horizontal  strip  of  short 
cells.)  The  size  of  these  square  cells  is  chosen  so  that  the  initial 
fireball  radius  is  at  least  8 cells  long.  The  time  at  which  the  fire- 
ball radius  will  expand  to  16  cells  is  then  estimated  by  means  of  an 
empirical  dimensionless  curve  of  fireball  radius  vs.  time.  At  that 
time  we  convert  all  variables  to  a new  grid,  having  cells  in  the  fire- 
ball region  that  are  twice  as  large  in  each  direction.  This  process 
is  repeated  as  needed. 

The  outer  part  of  the  grid  is  nonuniform,  each  cell  being 
20  per  cent  longer  or  wider  than  its  inner  neighbor.  Shock  waves  are 
allowed  to  run  off  the  edge  of  the  grid  when  they  become  moderate  to 
weak. 


The  artificial  diffusion  has  the  same  form  as  in  Ref.  1, 

Eq.  (42);  the  coefficients  qi , qz,  qs,  and  q^  are  chosen  to  be  0.5 
while  the  shock  is  strong  and  to  be  0.03  at  late  time.  The  artificial 
viscous  pressure  Q (Ref.  1,  Section  III)  is  unchanged  in  form;  its 
coefficient  qzz  is  chosen  to  be  1.0. 

D.  CALIBRATION  OF  RADIATION  LOSS 

Predicting  the  total  radiation  loss  for  low-altitude  bursts 
is  especially  important,  for  the  late- time  hydrodynamics  depends  on 
the  amount  of  energy  remaining  in  the  burst  region.  We  checked  the 
radiation  loss  predicted  by  our  method  (with  calibration  factor 
Cj^  = 1.0  in  Eq.  29)  against  that  predicted  by  the  RADFLO  code.  Using 
a one-dimensional  spatial  grid,  the  agreement  is  satisfactory.  (The 
power  vs.  time  curves  generally  peak  somewhat  too  early  and  too  high). 


I 


Using  the  two-dimensional  grid  required  for  interacting-fireball  calcula- 
tions results  in  radiation  loss  predictions  that  are  too  low.  The  coarser 
cells  necessary  in  a two-dimensional  calculation  seem  mainly  responsible 
for  this  error.  To  correct  for  this  we  increase  the  radiation  loss  rate 
(Eq.  29)  by  a calibration  factor  Cj^.  The  predicted  radiation  loss  for 
an  event  of  small  yield  is  satisfactory  if  Cj^  = 2,5  when  the  grid  is 
two-dimensional  and  c^^  = 1.0  when  it  is  one-dimensional. 
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IV  CONCLUSION 


The  code  has  been  used  to  calculate  a series  of  nine  problems 
involving  two  identical  bursts  separated  in  space  and  time.  The  results 
will  be  reported  elsewhere  by  Christian,  Stoeckly,  and  Schlueter. 


On  a GDC  7600  computer  using  the  RUN  compiler,  computation 
time  is  about  1 millisec  per  cell  per  time  step  for  hydrodynamics  and 
O.S  millisec  per  cell  per  time  step  for  radiation  diffusion.  After 
2 sec  of  problem  time  the  radiation  diffusion  terms  are  dropped  from 
the  calculations.  Total  computation  time  for  a problem  of  tv;o  close 
bursts  over  the  time  interval  0<t<120  sec  (2000  cycles)  using  a 
coarse  spatial  grid  averaging  20^40  cells  is  40  min  on  a CDC  7600. 
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